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We consider the problem of constructing optimal designs for pop- 
ulation pharmacokinetics which use random effect models. It is com- 
mon practice in the design of experiments in such studies to as- 
sume uncorrelated errors for each subject. In the present paper a 
new approach is introduced to determine efficient designs for nonlin- 
ear least squares estimation which addresses the problem of corre- 
lation between observations corresponding to the same subject. We 
use asymptotic arguments to derive optimal design densities, and the 
designs for finite sample sizes are constructed from the quantiles of 
the corresponding optimal distribution function. It is demonstrated 
that compared to the optimal exact designs, whose determination is 
a hard numerical problem, these designs are very efficient. Alterna- 
tively, the designs derived from asymptotic theory could be used as 
starting designs for the numerical computation of exact optimal de- 
signs. Several examples of linear and nonlinear models are presented 
in order to illustrate the methodology. In particular, it is demon- 
strated that naively chosen equally spaced designs may lead to less 
accurate estimation. 

1. Introduction. The work presented in this paper is motivated by some 
problems encountered in the design of experiments in a clinical trial to es- 
tablish the pharmacokinetics of Uzara® , a digitoxin related herbal diarrhea 
medication [based on Thiirmann, Neff and Fleisch (2004)]. These kinds of 
trials pose methodological design challenges because they require the esti- 
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mation of global population parameters in the presence of correlated mea- 
surement errors. The trial in question included a number of patients each 
given an oral application of Uzara® as well as, after a washout period, an 
intravenous application of digoxin (Lanicor®), where in both cases the re- 
sulting serum concentration of digitoxin was measured repeatedly during 
the next 36 hours. 

The relation between the time and the concentration in the analysis of the 
Uzara® trial can be described using the theory of one-compartment models 
with oral and, respectively, intravenous applications [Atkinson et al. (1993), 
Shargel (1993)]. In the intravenous case, the medication reaches the maxi- 
mum concentration in the blood almost immediately, and, after that, it is 
gradually eliminated from the body over time. Thus, the digitoxin concen- 
tration is modeled using the exponential elimination model 

(1.1) 7 1 (t,b) = he- b2t . 

In the case of the oral application, there is a gradual build up of the con- 
centration in the blood as the medication is absorbed through the digestive 
tract, while there is a simultaneous elimination process of the medication 
in the blood. Therefore, the concentration function is the solution of the 
differential equation of these two parallel processes. The resulting function 

(1.2) V (t,b) = b 3 (e- ht -e- b ^) 

is known as the (3 parameter) Bateman function [Garrett (1994)]. In both 
models, the function rj(t, 6) denotes the concentration function, t is the time 
(in hours), b= (61,62) and, respectively, b= (61,62,63) are the vectors of 
parameters. The parameters are assumed to vary between patients and the 
aim of the experiment is to estimate their global means (and sometimes 
variances) over all patients. 

Measurements within the same patient are usually correlated, and we 
assume this correlation to be proportional to the time lag between measure- 
ments, which is plausible as the random errors are usually caused by tempo- 
rary changes in the patients physical condition. Measurements for different 
patients are assumed to be independent. In the trial at question Thuermann 
considered n = 15 (oral), respectively, n = 14 (intravenous) measurements 
each on K = 18 patients. After a preliminary discussion with experts the 
measurements were taken at nonoptimized time points. An approximation 
of the covariance of a single patient can then be expressed as 

to, ^^m\^> +Vs , 

where rj(t,b) = (rj(ti, 6), . . . , rj(t n , b)) T denotes the vector of expected re- 
sponses at t\,...,t n and V e is the covariance matrix corresponding to this 
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data. This expression includes two sources of variation, the usual variation 
V e caused by random errors as well as the additional variation V p due to the 
random effect. 

Situations of this kind are rather common in the evaluation of the phar- 
macokinetics and the pharmacodynamics of drugs [see Buelga et al. (2005), 
Colombo et al. (2006), among others]. The corresponding processes are usu- 
ally modeled by linear or nonlinear random effects models, which try to esti- 
mate the population parameters, that is, the mean and the inter-individual 
variability of the parameters. Under the additional assumption of a normal 
distribution, the population characteristics are usually estimated by maxi- 
mum likelihood methods. In many cases, the likelihood cannot be evaluated 
explicitly and approximations are used to calculate the estimate. Efficient 
algorithms for estimation are available for this purpose [see Aarons (1999)]. 
Loosely speaking, under a Gaussian assumption this approach corresponds 
to weighted nonlinear least squares estimation. It was pointed out by sev- 
eral authors that the application of an appropriate design in these studies 
can substantially increase the accuracy of estimation of the population pa- 
rameters. Usually, the construction of a good design is based on the Fisher 
information matrix which cannot be derived explicitly in pharmacokinetic 
models with random effects. For this reason, many authors propose an ap- 
proximation of the likelihood [see, e.g., Mentre, Mallet and Baccar (1997), 
Retout, Mentre and Bruno (2002), Retout and Mentre (2003), Schmelter 
(2007a, 2007b), among others], which is used to derive an approximation 
for the Fisher information matrix. This matrix is considered in various op- 
timality criteria, which have been proposed for the construction of optimal 
designs for random effect regression models. 

In the present paper the investigation is motivated by the following is- 
sues. First, the estimation of the population mean and the construction 
of corresponding optimal designs for population pharmacokinetics strongly 
depends on the Gaussian assumption, which is usually made for compu- 
tational convenience. Moreover, the maximum likelihood estimates may be 
inconsistent if the basic distributional assumption is violated. As a con- 
sequence, the derived optimal designs might be inefficient. Second, most 
authors derive the approximation for the Fisher information matrix under 
the additional assumption that the random errors corresponding to the mea- 
surements of each individual are uncorrelated [see, e.g., Retout, Duffull and 
Mentre (2001), Retout, Mentre and Bruno (2002) and Retout and Mentre 
(2003), among many others]. However, this assumption is not realistic in 
many applications of population pharmacokinetics. Thus, a general concept 
for constructing optimal designs in the presence of correlated observations 
is still missing. Third, even if the Gaussian assumption and the assumption 
of uncorrelated errors for each subject can be justified, the numerical con- 
struction of the estimate and the corresponding optimal design is extremely 
hard. 
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In the present paper we address these issues. To derive a good design, 
we consider nonlinear least squares estimation in random effect regression 
models. Note that this estimation does not require a specification of the 
underlying distributions. For this estimation, we introduce a methodology 
which can be used to derive efficient or optimal designs in very general 
situations. More precisely, we embed the discrete optimal design problem in 
a continuous optimal design problem, where a nonlinear functional of the 
design density has to be minimized or maximized. This approach takes into 
account the correlation dependence and yields an asymptotic optimal design 
density, which has to be determined numerically in all cases of practical 
interest. For finding the optimal density, we propose an algorithm based on 
polynomial approximation. For a fixed sample size and for each individual, 
an exact design can be obtained from the quantiles of the corresponding 
optimal distribution function. It is demonstrated by examples that these 
designs are very efficient. Moreover, the designs, derived from the asymptotic 
optimal design density, are very good starting designs for any procedure of 
local optimization for finding the exact optimal designs. To our knowledge, 
the proposed method is the first systematic approach to determine optimal 
designs for linear and nonlinear mixed effect models with correlated errors. 

The remaining part of this paper is organized as follows. In Section 2 we 
consider the case of a linear random effect model and explain the basic design 
concepts. In Section 3 we introduce the approach for deriving the asymptotic 
optimal designs for linear regression models with correlated observations by 
employing results of Bickel and Herzberg (1979). In Section 4 we present 
results for nonlinear regression models with correlated random errors and 
derive D-optimal designs and optimal designs for estimating the area under 
the curve in the compartmental model. Finally, the Uzara® and Lanicor® 
trials are re-analyzed and optimal designs for the model (1.2) are determined. 
In particular, we show that the design proposed by the experts is rather 
efficient, while a naively chosen equidistant design can yield substantially 
less accurate estimates. 

2. Statement of the problem. Consider the common random effect linear 
regression model 

(2.1) Y ij = bjg(t ij ) + e ij , i = 1, . . . ,K; j = 1, . . . ,m, 

where Yij denotes the jth observation of the ith subject at the experimen- 
tal condition tij, ■ ■ ■ ,£R,n K are centered random variables with vari- 
ances depending on t, Var(ejj) = a 2 h 2 (tij) for some positive function h(t), 
g(t) = (gi(t), . . . ,g p (t)) T is a given vector of linearly independent regression 
functions, and hi is a p-dimensional random vector representing the individ- 
ual parameters of the ith subject, i = 1, . . . ,K . The explanatory variables 
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can be chosen by the experimenter from a compact interval T. We assume 
that errors £j = (en, . . . , £j n J for different subjects are independent, but the 
errors for the same subject are correlated, that is, 

(2.2) Cov(eij,s is ) = cr 2 /i(%)/i(t is )(7r(% - t is ) + (1 - j)S jtS ), 

where 7 £ [0,1] is a constant, r(t) is a given correlation function such that 
r(0) = 1, and 5j s denotes Kronecker's symbol. Let V e be the corresponding 
covariance matrix. Assume that the individual parameters 6j are drawn from 
a population distribution with mean (3 and covariance matrix V p , and they 
are independent of the random variables £j. This means that the covariance 
between two observations at the time tij and the time ti s (j 7^ s) is 

Cov(Yij,Y is ) = g T (tij)Vpg(t is ) + a 2 h(tij)h(t is )^r(tij - t is ), 

while the variance of Yij is given by g 1 ' (tij)V p g(tij) + a 2 h 2 (tij) . It was shown 
by Schmelter (2007b) that an optimal design necessarily advises the exper- 
imenter to perform observations of all subjects at the same experimental 
settings, that is, Uj = tj (i = 1, . . . ,K, j = 1, . . . , n). Consequently, we define 
an exact design £ = {t±, . . . ,t n } as an n-dimensional vector which describes 
the experimental conditions for each subject. Without loss of generality, we 
assume that the design points are ordered, t\ < ■ ■ ■ <t n . 

Suppose that n observations are taken according to the design £. Then 
the model (2.1) for the ith subject can be written as 

(2.3) Y i = X g b i + e i ; i = l,...,K, 

where Y i = (Y a ,..., Y in ) T , the matrix X g is given by X g = (#(ii), . . . , g{t n )) T , 
and E; is now a centered random variable with variance a h 2 (U). This model 
is a special case of the random-effect models discussed in Harville (1976), 
which are called generalized MANOVA. According to Fedorov and Hackl 
(1997), the (ordinary) least square estimate of (3 minimizes 

K n 

1=1 j=i 

Define f(t) = g(t)/h(t) and X = (/(ti), . . . , f(t n )) T . Then the covariance 
matrix of the ordinary least squares estimate /3ols is given by 

D(/3 OLS ) = ^(x T xy 1 x T (v £ + XV P X T )X(X T X)- 1 

(2.4) K 

= -((x T xy 1 x T v £ x(x T x)- 1 + V p ). 

Alternatively, if the covariance matrix V e of the errors and the covariance 
matrix of the random effects V p were known (or can be well estimated), the 
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(weighted) least squares statistic 
1 K 

(2.5) /3 W ls = ^ ]T(X T (F e + XV p X T )- l X)~ l X T (V e + XVpX^Yi, 

i=i 

can be used to estimate the parameter (3. The covariance matrix of the 
estimate /3\vls is given by 

(2.6) D(/3 W ls) = ^(X T (V £ +XV p X T r 1 X)-\ 

Since the expression (2.4) is simpler than (2.6) (which requires two dif- 
ferent inversions), the design methodology developed in this paper is based 
on the covariance matrix of the ordinary least squares estimate. Addition- 
ally, we will demonstrate that the optimal designs obtained by minimizing 
functionals of the covariance matrix of the ordinary least squares estimate 
are also very efficient for weighted least squares estimation. 

We call a design optimal design if it minimizes an appropriate functional 
of the covariance matrix of the least squares estimate. Since we consider 
optimality criteria that are linear with respect to scalar multiplication of 
the covariance matrix, we put K = 1 without loss of generality. 

3. Asymptotic optimal designs. Although the theory of optimal design 
has been discussed intensively for uncorrelated observations [see, e.g., Fe- 
dorov (1972), Pazman (1986) and Atkinson and Donev (1992)], less results 
can be found for dependent observations. For linear and nonlinear random 
effect models, optimal designs under the assumption of uncorrelated errors 
have been investigated in Schmelter (2007a, 2007b), Mentre, Mallet and 
Baccar (1997) and Retout and Mentre (2003), among others. For fixed effect 
regression models with the presence of correlated errors, it was suggested to 
derive optimal designs by asymptotic considerations. Sacks and Ylvisaker 
(1966, 1968) have considered a fixed design space, where the number of 
design points in this set tends to infinity. As a result, the asymptotic op- 
timal designs depend only on the behavior of the correlation function in a 
neighborhood of the point 0. In the present paper we use the approach of 
Bickel and Herzberg (1979) and Bickel, Herzberg and Schilling (1981), who 
have considered a design interval expanding proportionally to the number 
of observation points. This case is equivalent to the consideration of a fixed 
interval with the correlation function depending on the sample size. To be 
precise, we assume that the design space is given by an interval T and the 
design points of a sequence of designs £ n = {t\ n , . . . , t nn } are generated by a 
function a(-) in the form 



(3.1) t jn = a((j-l)/(n-l)), j = l,...,n; 
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and a : [0, 1] — > T denotes the inverse of a distribution function. Note that the 
function a(-) is obtained from the density of the weak limit of the sequence 
£ n as n — > oo. For example, if T = [—1, 1], the function a(u) = (2u — 1) cor- 
responds to the equally-spaced designs with distribution function a~ 1 (x) = 
and density (a~ 1 )'(x) = ^I\-i t i](x). Furthermore, we assume that the 
correlation function r(t) of the errors ej in (2.2) depends on n in the form 

(3.2) r n (t) = p{nt), 

such that p(t) = o(t) as t — > oo. For the numerical construction of asymptotic 
optimal designs, we derive an asymptotic representation for the covariance 
matrix of the least squares estimate in the following lemma. For this purpose, 
we impose the following regularity assumptions: 

(CI) The functions fi(t), . . . , f p (t) are linearly independent, bounded on the 
interval T and satisfy the first order Lipschitz condition: 

\fi(t)-fi(s)\<M\t-s\ and 

\fi{t)\<M for &l\t,sET,i = l,...,p. 

(C2) The function a(-) is twice differentiable and there exists a positive 
constant M < oo such that for all u G (0, 1), 

(3.3) -7j<a'(«)<M, \a"(u)\<M. 

(C3) The correlation function p is differentiable with bounded derivative 
and satisfies p'{t) < for sufficiently large t. 

Assumption (CI) refers to the continuity of the response as a function of t 
and it is satisfied for most of the commonly used regression models. More- 
over, most of the commonly used correlation structures satisfy assumption 
(C3) on a compact interval. Finally, assumption (C2) refers to the charac- 
teristics of a design, requiring, loosely speaking, that observations should 
not be clustered. This is quite reasonable due to ethical aspects. The follow- 
ing result is obtained by similar arguments as given in Bickel and Herzberg 
(1979) and, hence, its proof is omitted. 

Lemma. Assuming that conditions (CI), (C2) and (C3) are satisfied, 
then the covariance matrix of the least squares estimate have the form 

(3.4) D(/3 OLS ) = —(W-\a) + 2 1 W- l (a)R(a)W-\a)) + V p + o(-), 

n \n J 

where the matrices W(a) and R(a) are defined by 

W(a)=( f 1 / i (a(«))/,(a(«))d«V , 
\Jo J i,j=l 

R(a) = ( [ /i(a(u))/ i (a(«))Q(a'(«))^y 

V0 / i,j=i 
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and the function Q(-) is given by Q(t) = Yl^=i P(jt)- 

Note that only the first term in (2.4) [and (3.4)] depends on a (asymptotic) 
design, and we propose to use this term for the construction of optimal 
designs. If the function a(-) is the inverse of a continuous distribution with 
density, say, tp, then a'(t) = l/p(t) and, for large n, the first term of the 
covariance matrix can be approximated by the matrix V(ip)/n, where the 
p x p matrix V is given by 

(3.5) V(ip) := a 2 {W- x {<p) + 2 1 W- 1 {p)R{p)W- 1 {p)). 
In (3.5) the matrices W(<p) and R(<p) have the form 

w(<p) = ( I fiWsitMQdt) , 

RifP) = ( ! m)f j {t)Q{l/i P {t))y{t)dt) P . 

\JT J i,j=l 

A density will be called an asymptotic optimal density ip* , if it minimizes an 
appropriate functional of the matrix V((p). Numerous criteria have been pro- 
posed in [Silvey (1980), Atkinson and Donev (1992), Pukelsheim (1993)] and, 
exemplarily, we consider the D- and c-optimality criteria which minimize 
— det V(ip) and c T V(ip)c for a given vector c € W 3 , respectively. The applica- 
tion of the proposed methodology to other optimality criteria is straightfor- 
ward. The general procedure for constructing an efficient design minimizing 
a given functional of the covariance matrix of the least squares estimate is 
as follows: 

(1) Specify the correlation function p(-) in (3.2) and compute Q(-). 

(2) Compute the asymptotic optimal design density <p* that minimizes an 
appropriate functional of the matrix V(ip) in (3.5). 

(3) Derive the exact design for a fixed sample size n by calculation of the 
quantiles of the distribution function $* that corresponds to p* , namely, 

(3.6) t . n = ($*)-! ^izlY i = l,..., n . 

The optimal density ip* in step (2) can be determined as follows. We consider 
the parametric representation of a density by a polynomial in the form 

,q = (Po+Pi^H h£vQ+ 

j T (Po +Pit H \-p r t r )+dt 

and apply the Nelder-Mead algorithm to find the optimal density that min- 
imizes the specified functional of the matrix V(p) with respect to po, . . . ,p r . 
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One can run the algorithm for different degrees of the polynomial and differ- 
ent initial values and choose the density corresponding to the minimal value 
of the optimality criterion. All integrals can be calculated by the Simpson 
quadrature formula. We found that the minimal value of the criterion is 
negligibly decreasing for polynomials of degree larger than r = 6. We also 
investigated the case where a density is represented in terms of rational, ex- 
ponential or spline functions. The results were very similar and, on the basis 
of our numerical experiments, we can conclude that the optimal density ip* 
can be very well approximated by 6-degree polynomials. 

The derived designs from the asymptotic theory can be used to construct 
efficient designs for a given sample size as specified in step (3) of the algo- 
rithm. Alternatively, exact optimal designs can be determined by employing 
the derived designs as initial values in a (discrete) optimization procedure. 
More precisely, for the determination of an exact optimal design, the above 
procedure can be extended by the fourth step: 

(4) Determine an exact optimal design, that minimizes a functional of the 
covariance matrix in (2.4), by using the Nelder-Mead algorithm with an 
initial n-point design, which has been found in step (3). 

We illustrate this methodology for a quadratic regression model by the fol- 
lowing example. Also, in Section 4 we extend the approach for designing 
for nonlinear models and investigate its performance for the compartmen- 
tal model. In particular, we demonstrate that the designs derived from the 
asymptotic consideration are very efficient compared to exact optimal de- 
signs. 

Example 1 (Quadratic regression). We illustrate the proposed approach 
for constructing optimal designs for the classical quadratic regression model 
with homoscedastic errors. For this model, we have p = 3, fi(t) = 1, fi(f) = t, 
/2W = t 2 and T = [—1, 1]. Let the correlation function in (3.2) be given by 
p(t) = e~ xt and 

(3.7) r n (t) = e- Xnt , A > 0. 

The asymptotic D-optimal densities for different choices of the parameters 
A and 7 are shown in Figure 1. Note that the numerically calculated optimal 
design densities are symmetrical, but we were not able to prove the symmetry 
of the asymptotic optimal density. We observe that the L)-optimal density 
converges to the density of the uniform design, if 7 — > 1 or A — > 0. On the 
other hand, if 7 — > or A — > 00, it can be seen that the asymptotic D-optimal 
design density is more concentrated at the points —1, and 1, which are 
the points of the exact D-optimal design for the quadratic fixed effect model 
with uncorrelated observations [see Gaffke and Krafft (1982)]. Such behavior 
is natural because the errors are less correlated as 7 — > or A — > 00. 
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Fig. 1. Asymptotic D-optimal design densities for least squares estimation in a random 
effect quadratic model for different choices of parameters in the covariance function (2.2) 
with r(t) defined by (3.7). Right panel 7 = 0.6; left panel A = 0.2. 



Further, we investigate the efficiency of an exact design derived from the 
asymptotic theory for least squares estimation, where the parameters in the 
correlation function (2.2) are given by 

(3.8) 7 = 0.6, a 2 = 0.5, F p = diag(a^,^ 2 ) = diag(0.3 2 ,0.3 2 ). 

Let £JJ be an n-point equidistant design and £° be an n-point design obtained 
by the transformation (3.6) from the asymptotic optimal density. The points 
of the design and the exact optimal designs for least squares estimation 
for A = 1.2 are displayed in the left part of Figure 2. In the right part of 
Figure 2, we present the efficiencies for ordinary least squares estimation, 

det[(X^)-i X T FeX(x T x) -i + v] x VP 

effoLs(f) " 



.det[(x5 LS XoLs)- 1 ^SLS^^OLs(^LS^OLs)- 1 + V p ] 
and for weighted least squares estimation, 

/ det[X T (V £ + XV p X T )^X] \^ 

WLSlU UetiXT LS (V £ + X WLS V p XT LS )-iX WLS ]) 

for two designs: the design derived from asymptotic theory, and the 
uniform design Here X denotes the design matrix obtained for the design 
£ under investigation, while Xols an d -A^wls correspond to the optimal 
exact design for ordinary and weighted least squares estimation respectively. 
We observe that the design points concentrate in three regions containing 
the points of the exact D-optimal design for a quadratic regression with 
uncorrelated errors. It is also noteworthy that the designs derived from the 
asymptotic theory are very efficient, in particular, for weighted least squares 
estimation. 
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Fig. 2. Left part: Various designs for ordinary and weighted least squares estimation 
in the random effect quadratic regression model. Exact D-optimal designs derived from 
asymptotic theory: ball; exact D-optimal for ordinary least squares estimation: diamond; 
exact D -optimal designs for weighted least squares estimation: triangle. Right part: Effi- 
ciency of the designs and equidistant design for ordinary and weighted least squares 
estimation. The parameters are given by (3.8), where A= 1.2. 

4. Nonlinear random effect models. In this section we extend the method- 
ology to the case of nonlinear random effect models, which have found con- 
siderable interest in the literature on pharmacokinetics. In this case, the 
results of experiments are modeled by 

(4.1) Yi j = r](tj,bi) + £ ij , i = l,...,K;j = l,...,n. 

Since the model (4.1) is nonlinear with respect to the variables 6j, there is 
no analytical expression for the likelihood function and various approxima- 
tions have been considered in the literature [see Mentre, Mallet and Baccar 
(1997), Retout, Mentre and Bruno (2002), Retout and Mentre (2003), among 
others] . These approximations are used for the calculation of the maximum 
likelihood estimate and the corresponding Fisher information matrix. Al- 
ternatively, an estimate of the population mean f3 can be obtained as an 
average of the nonlinear least squares estimates hi for the different individ- 
uals, but due to the nonlinearity of the model, an explicit representation of 
the corresponding covariance matrix cannot be derived. Following Retout 
and Mentre (2003), we employ a first-order Taylor expansion to derive an 
approximation of this covariance matrix. To be precise, we obtain (under 
suitable assumptions of differentiability of the regression function) the ap- 
proximation 

(4.2) r,(t,b)^ V (t,f3)+g(t,f3)(b-i3) T , 



where 



12 



H. DETTE, A. PEPELYSHEV AND T. HOLLAND-LETZ 



denotes the gradient of the regression function with respect to b. This means 
that the nonlinear model (4.1) is approximated by the linear model (4.2). 
For the construction of the optimal design, we assume that knowledge about 
the parameter ft is available from previous or similar experiments. This cor- 
responds to the concept of locally optimal designs, introduced by Chernoff 
(1953) in the context of fixed effect nonlinear regression models. Usually, 
locally optimal designs serve as benchmarks for commonly used designs and 
are the basis for the construction of optimal designs with respect to more 
sophisticated optimality criteria including the Bayesian and minimax ap- 
proach [see Chaloner and Verdinelli (1995) or Dette (1995)]. 

Following the discussion in Section 3, we define the function f(t,b) = 
g(t, b) /h(t) to account for heteroscedasticity. Note that the covariance matrix 
of the nonlinear least squares estimate in the model (4.1) is approximated 
by replacing the matrix X in model (2.1) with f{t) = f(t, b)\b=/3, and the 
methodology described in Sections 2 and 3 can be applied to determine 
efficient designs. In the context of dose finding studies, it has been shown 
by means of a simulation study that the approximation (4.2) has sufficient 
accuracy for the construction of optimal designs [see Section 5 in Dette et al. 
(2008) for more details]. 

We further illustrate this concept by giving several examples for the case of 
homoscedastic errors. First, we investigate D- and c-optimal designs for the 
random effect model, which has recently been studied by Atkinson (2008). 
Next, we re-analyze the Uzara® and Lanicor® trials introduced in Section 1. 

4.1. D-optimal design for a random effect compartmental model. We 
consider the random effect compartmental model with first-order absorp- 
tion, 

(4.3) Tf (t,b) = -\-(e- bat -e-^ t ). 

bi - b 2 

The model (4.3) is a special case of the Bateman function, defined in the 
introduction [see Garrett (1994)], and has found considerable attention in 
chemical sciences, toxicology and pharmacokinetics [see, e.g., Gibaldi and 
Perrier (1982)]. The optimal design problem in the compartmental fixed ef- 
fect model has been studied by numerous authors [see, e.g., Box and Lucas 
(1959), Atkinson et al. (1993), Dette and O'Brien (1999), Biedermann, Dette 
and Pepelyshev (2004), among others], but much fewer results are available 
under the assumption of random effects. Recently optimal approximate de- 
signs for a random-effect compartmental model (4.3) have been determined 
by Atkinson (2008), but we did not find results about exact designs for these 
models in the presence of correlated errors. In the present paper we derive 
such designs from the asymptotic optimal design density and compare these 
designs with the exact optimal designs. 
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Fig. 3. Asymptotic D-optimal design densities for nonlinear least squares estimation in 
the compartmental model (4.3) for different choices of the parameters in the covariance 
function (2.2) with r(t) defined by (3.7). Left part: A = 0.2; right part: 7 = 0.6. 



Note that the gradient of the function rj(t, b) with respect to b is given by 
, ,x fb 2 (e- bix -e- b2X ) + (b 1 2 x-b 1 b 2 x)e- bia 

9{tM = { (p^w 

(4.4) 

fri(e- bix - e~ b2X ) + (bi 2 x - hb 2 x)e- b2 

(h-b 2 y 

In order to illustrate the methodology, we consider the same scenario as in 
Atkinson (2008) and assume that the parameters of the population and error 
distributions are the following: = (1,0. 5) T , 

(4.5) 7 = 0.6, a 2 = 0.01, V p = diag(a| 1 , a%) = diag(0.1 2 , 0.05 2 ), 

and the design space is given by the interval T = [0, 10] . We assume that the 
function r(t) in (2.2) is given by (3.7). The asymptotic D-optimal design 
densities for different values of the parameters are shown in Figure 3. We 
observe that, for 7 — > 1 or A — > 0, the D-optimal design densities approximate 
the uniform design, while, for larger values of A or small values of 7, the 
asymptotic D-optimal designs put more weight at two specific regions of the 
design space. This fact corresponds to intuition, because the (approximate) 
D-optimal design for the model (4.3) with uncorrelated observations is a 
two-point design [see, e.g., Box and Lucas (1959)]. 

Further, we investigate the performance of the uniform and an exact de- 
sign derived from asymptotic theory. For this purpose we define £JJ as an 
n-point equidistant design {10/n, 20/n, . . . , 10} and £JJ as the ra-point design 
obtained by the transformation 




t j = ($T 1 (3/(n+l)), j = l,-..,n, 
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where $* denotes the distribution function corresponding to the asymptotic 
D-optimal design density. Note this transformation is slightly different from 
the transformation (3.1) in order to exclude the point from the design 
points. Obviously, it is not reasonable to take observations t = in model 
(4.3), because it is assumed that the drug is administered at time t = 0. 
The corresponding points of the exact designs are depicted in the left part 
of Figure 4, while the right part of the figure shows the D-efficiencies of the 
different designs. We observe that the designs derived from the asymptotic 
theory have a substantially larger D-efnciency compared to the uniform 
design. For example, for n = 6, the D-efficiency of the uniform design is 
approximately 50% for ordinary and weighted nonlinear least squares esti- 
mation, while the Z)-efficiency of the design £° is close to 90%. 

It is worthwhile to mention that, in nonlinear random effect models, the 
optimal designs depend additionally on the mean j3 of the distribution of the 
population parameters 6j. Therefore, it is of interest to investigate the sensi- 
tivity of the designs with respect to a misspecification of this parameter. For 
a study of the impact of such a misspecification on the efficiency of the result- 
ing designs, we consider the case n = 4 and n = 6 and the corresponding de- 
signs ££ = {1.04, 2.01, 3.16, 4.33} and = {0.83, 1.47,2.32,3.30,4.20,5.20}, 
respectively. In Figure 5 we display the efficiencies 

-i/p 



(4.6) 



det[(X T X)~ 1 X T V e X(X T X)- 1 + V p ] 



det [( X OLS,/3 X OLS,/3) 1 Xg LS/3 V r £ X LS,/3(^OLS,/3 X OLS,/3) 1 + Vp. 



for different values of /3. Here X denotes the design matrix obtained from 
the design under the assumption that = (1,0, 5) T , while the matrix 




5 6 7 



Fig. 4. Left part: Various designs for ordinary and weighted nonlinear least squares es- 
timation in the compartmental model (4.3). Exact D-optimal designs derived from asymp- 
totic theory: ball; exact D-optimal for ordinary least squares estimation: diamond; exact 
D-optimal designs for weighted least squares estimation: triangle. Right part: Efficiency of 
the designs and for ordinary and weighted least squares estimation. The parameters 
are given by (4.5), where A = 0.2. 
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Fig. 5. D -efficiencies of the designs £4 = {1.04,2.01,3.16,4.33} (left part) and the design 
fg = {0.83, 1.47, 2.32, 3.30, 4.20, 5.20} (right part) (these designs are obtained from asymp- 
totic density) for nonlinear least squares estimation in the random effect compartmental 
model (4.3) , if the mean of the population distribution has been misspecified. The param- 
eters of the population distribution are given by (4.5) with A = 0.2. 



Xols,/3 corresponds to the exact D-optimal design for nonlinear least squares 
estimation for a specific (3. The efficiencies are plotted in Figure 5 for the 
rectangle [fif ) - 3a ( ) , pf ) + 3<r ( 0) ] = [0.7,1.3] x [0.35,0.65]. It can be seen 
that the exact optimal design, derived from the asymptotic theory, has a 
high L>-efficiency with respect to misspecification of the population mean /? 
over a broad range. 

4.2. Optimal designs for estimating the AUC. In some bioavailability 
studies, the aim of experiments is the estimation of the area under curve 

AUC= / rj(t,P)dt. 
Jo 

For the compartmental model (4.3), we obtain AUC = 1/&2- It can be shown 
that the locally AUC-optimal design for model (4.3) minimizes the variance 
of the nonlinear least squares estimate for the parameter 02- This variance 
is approximately proportional to 

(o, i)([x T x)- 1 x T v £ x{x T xy l + v p )(o, if. 

This expression corresponds to the c-optimality criterion, which has been 
discussed extensively in the literature for fixed effect models with uncorre- 
cted observations [see, e.g., Ford, Torsney and Wu (1992), Fan and Chaloner 
(2003) and Dette et al. (2008), among others]. The asymptotic optimal de- 
sign densities for estimating the area under the curve are shown in Figure 
6. We observe again that the optimal density is close to the uniform design 
density if A — >■ or 7 — > 1. On the other hand, if A is large or 7 — > 0, the 
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Fig. 6. Asymptotic optimal densities for estimating the area under the curve in the com- 
partmental model (4.3) for different choices of the parameters in the covariance function 
(2.2) with r(t) defined by (3.7). Left part: A = 0.2, right part: 7 = 0.6. 

AUC-optimal design density has a narrow support. This fact reflects that 
the optimal design for estimating the area under the curve in the fixed effect 
compartmental model with uncorrelated observations is a one-point design. 
In Figure 7 we show the designs derived from the asymptotic optimal design 
densities, and the exact optimal designs for estimating the area under the 
curve in the compartmental model. We observe that the designs derived 
from the asymptotic optimal design density, are very close to the exact 
optimal designs for least squares estimation of the area under the curve. 
Moreover, the design yields a substantial improvement in efficiency com- 
pared to the uniform design. Again we observe that the designs derived for 
ordinary least squares estimation also have excellent efficiencies for weighted 
least squares estimation of the AUC (see the right part of Figure 7) . 

4.3. Optimal designs for estimating the A UC in the Uzara® and Lanicor® 
trials. In this section we consider the optimal design problem for estimat- 
ing the AUC in the two examples presented in the introduction. The med- 
ical background of these examples was a small pilot trial of 4 patients 
[Thiirmann, Neff and Fleisch (2004)], where it was observed that patients 
taking the over-the-counter herbal diarrhea medication Uzara® (in the form 
of drops, i.e., oral application) showed high values in medical assays designed 
to measure the blood serum concentration of digitoxin, a potent treatment 
against heart insufficiency. This is caused by the chemical similarity of these 
two substances and can result in major complications in establishing treat- 
ment programs for heart insufficiency. It was thus decided to compare the 
pharmacokinetic properties of an oral application of Uzara® to the proper- 
ties of the usual intravenous application of a regular digitoxin medication 



OPTIMAL DESIGNS FOR RANDOM EFFECT MODELS 



17 




Fig. 7. Left part: Various designs for estimating the area under the curve in the compart- 
rnental model (4.3). Designs derived from asymptotic theory: ball; exact optimal designs 
for least squares estimation of the area under the curve: diamond; exact optimal designs 
for weighted least squares estimation: triangle. Right part: Efficiency of the designs £JJ and 
£,n f or ordinary and weighted least squares estimation of the AUG The parameters are 
given by (4.5), where A = 1.2. 



(Lanicor®) on a larger sample size of 18 patients, with 15 measurements 
each. The main focus of the comparison was the area under the concentra- 
tion curve as a measure of the total effect of an application. A preliminary 
design was proposed by experts in order to allow precise estimation of this 
property, and the study was carried out according to this design. We will 
now investigate the efficiency of this and a naively chosen equally spaced 
design compared to a design generated using the methods presented in this 
paper. 

In order to compute the design densities, it is necessary to estimate the 
parameters for both of the models. To do so, we used the full 18 patient data 
set, however, this could have been done as well using only the 4 patients of 
the pilot study. We estimated parameters using a combination of maximum 
likelihood and least squares techniques [see Pinheiro and Bates (1995)]. For 
the intravenous part corresponding to model (1.1), we have received the 
estimates /3 = (30, 0.75) T , 

y - ( 9 °- 189 ^ 
P V 0.189 0.0049 7 

for the population parameters, while the estimates of the parameters in the 
covariance matrix (1.3) are given by 7 = 0.8, A = 0.05, a 2 = 0.6. 
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For the oral application corresponding to model (1.2), the estimates of 
the parameters are given by (3 = (0.2, 0.135, 28) T , 



Vn 




0.0019 
0.0016 




and 7 = 0.8, A = 0.01, a = 0.2. Note that the entries in the positions (1, 3), 
(2,3), (3,1) and (3,2) are not identical 0, but smaller than 10~ 5 . Based 
on a preliminary discussion with experts, the physicians decided to take 
measurements for both experiments at nearly identical (nonoptimized) time 
points 



and 



0, 0.5, 1, 1.5, 2, 3, 4, 5, 6, 8, 10, 12, 15, 24, 36 



0.25, 0.75, 1.25, 2, 3, 4, 5, 6, 8, 10, 12, 15, 24, 36, 



respectively. Note that in the second design an additional measurement was 
taken at t = —0.25, that is, before the intravenous injection. Since this point 
is out of the scope of the exponential evasion model, the observation at this 
time has been excluded from further considerations. Thus, the Lanicor® 
trial has n = 14 observations. 

For the estimated parameters, we have derived the asymptotic optimal 
design densities, which are depicted in the left and right part of Figure 8 for 
the Uzara® and Lanicor® trials, respectively, where the design interval is 
given by T = [0, 36]. The resulting designs from these densities are given by 

2.09, 4.55, 7.49, 10.8, 13.9, 16.8, 19.2, 21.5, 23.6, 25.7, 27.7, 29.8, 31.9, 34.0, 36 




15 20 25 30 



Fig. 8. Asymptotic optimal density for estimating the area under the curve in the 
^-parameter Bateman (left part) and exponential evasion model (right part). 
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for the Uzara® trial and 

0.432, 0.85, 1.25, 1.66, 2.06, 2.46, 2.87, 3.28, 3.7, 4.12, 4.55, 5.01, 5.54, 36 

for the Lanicor® trial. This means that the asymptotic optimal design for 
the oral application is close to an equally spaced design, while the optimal de- 
sign for the intravenous application is much more focused on measurements 
at early time points. This result is plausible in an exponential elimination 
model. Calculating the efficiencies for ordinary least squares estimation of 
these asymptotic optimal designs compared to the exact optimal design, we 
obtain efficiencies 0.99 and 0.95, respectively, indicating that these designs 
are quite good. 

Upon numerical studies, we can further conclude that the constructed 
designs are robust to misspecification of the initial guess of parameters. For 
example, variations of A (or any of the other parameters) by up to 50% yield 
a drop of efficiency by less than 0.02. 

Note that in the examples presented in Sections 3, 4.1 and 4.2, the ef- 
ficiencies of the derived designs for weighted least squares estimation are 
very similar to the efficiencies for ordinary least squares estimation, and 
a similar observation has been made for the two trials under investigation. 
The efficiencies for weighted least squares estimation of the designs based on 
the asymptotic density are 99% (Uzara®) and 97% (Lanicor®), even better 
than the efficiencies for ordinary least squares estimation. 

We now investigate the performance of the designs which were actually 
used in the clinical trial. We found that these designs have efficiency 0.96 
and 0.92 for estimating the AUC through OLS estimation in the Uzara® and 
Lanicor® trials, respectively. Thus, the preliminary designs, recommended 
by experts, are rather efficient in both trials. 

Let us now investigate the performance of naively chosen equidistant de- 
signs: the design 

0, 2.6, 5.1, 7.7, 10.3, 12.8, 15.4, 18.0, 20.6, 23.1, 25.7, 28.3, 30.8, 33.4, 36 

for the Uzara® trial, and the design 

0, 2.8, 5.5, 8.3, 11.1, 13.9, 16.6, 19.4, 22.2, 24.9, 27.7, 30.5, 33.2, 36 

for the Lanicor® trial. Comparing these designs to the optimal designs, we 
obtain efficiencies of 0.97 and 0.41 for ordinary least squares estimation of 
the AUC in the Uzara® and Lanicor® trials, respectively. 

It was pointed out by a referee that it is of some interest to investigate the 
performance of the optimal designs proposed in this paper for the estimation 
of the variance of the random effects (i.e., the parameters of the matrix V p ). 
These parameters are usually estimated by maximum likelihood techniques 
[see Retout and Mentre (2003), among others] and the corresponding infor- 
mation matrix of these estimates is of a fundamentally different structure 
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compared to the variance matrix of the least squares estimate. For the two 
optimal designs we have calculated the D-efficiencies for estimating the di- 
agonal elements of the matrix V p , which are 96% and 79% in the Uzara® 
and Lanicor® trial, respectively. The designs actually used in the trial have 
efficiency 98% and 85%, while the corresponding efficiencies of the uniform 
designs are 98% and 63%, respectively. Thus, the proposed optimal designs 
for AUC-estimation also have reasonable efficiency for estimation of the co- 
variance matrix of the population distribution. 

Summarizing the discussion in this example, we conclude that the equally 
spaced design in the Uzara® trial is very close to the optimal design deter- 
mined by the proposed methodology, and it is for this reason very efficient. 
However, the equally-spaced designs do not always have high efficiency. In 
the Lanicor® trial, the use of naively chosen designs yields considerably less 
accurate estimates. For this reason, the application of experimental design 
techniques in the context of pharmacokinetics trials is strictly recommended. 
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